
# map plot fun
map.plot.fun <- function(area.name){
  
  df3.shapes.temp <- df3.shapes[grep(area.name, df3.shapes$lsoa11nm),]
  years <- unlist(strsplit(df3.shapes.temp$years, ", "))

  plot(df3.shapes.temp[,"lsoa11cd"], 
       col = heat.colors(n.color.bins, rev = TRUE)[df3.shapes.temp$color.bins],
       main = paste0("Hate crimes per capita in ", area.name, ", ", min(years), " - ", max(years)), 
       axes = FALSE )
  legend("bottomright", cex = 0.8,
         title = "# crimes per 100k",
         fill = c(heat.colors(n.color.bins, rev = TRUE), "white"),
         legend = c(levels(df3.shapes$hcrime.bins), "Not observed"))
}

map.plotter <- function(area){
  for(i in 1:length(area)){
    main.lad.code <- names(table(df2.shapes$ladcode[grepl(area[i], df2.shapes$lsoa11nm)])[which.max(table(df2.shapes$ladcode[grepl(area[i], df2.shapes$lsoa11nm)]))])
    pdf(paste0("plots/", area[i], ".plot.pdf"), width = 8, height = 8)
    our.lad.plot.fun(area=area[i], main.lad.code=main.lad.code)
    dev.off()
  }
}


# prediction bootstrap:
boot.fun.internal <- function(the.models, x.levels, int){
n.xlevels <- nrow(x.levels)
pred.data <- data.frame(x.levels, 
                        fit.inc = NA, lwr.inc = NA, upr.inc = NA,
                        fit.env = NA, lwr.env = NA, upr.env = NA,
                        fit.edu = NA, lwr.edu = NA, upr.edu = NA)
data.row <- sample(1:nrow(df5), 1)

if(int == FALSE){
  new.data1 <- suppressWarnings(data.frame(x.levels, 
                                           lapply(the.models, function(x){df5[data.row, attr(x$terms, "term.labels")[2]]} ),
                                           df5[data.row, attr(the.models[[1]]$terms, "term.labels")[-c(1,2)]]))
  
  colnames(new.data1) <- c(c( do.call("rbind", lapply(the.models, function(x){attr(x$terms, "term.labels")[1:2]} )) ), attr(the.models[[1]]$terms, "term.labels")[-c(1:2)] )
  
  pred.data[,grepl("\\.inc", colnames(pred.data))] <- predict(the.models[[1]], newdata = new.data1, type = "response", interval = "confidence", level = 0.95)
  pred.data[,grepl("\\.env", colnames(pred.data))] <- predict(the.models[[2]], newdata = new.data1, type = "response", interval = "confidence", level = 0.95)
  pred.data[,grepl("\\.edu", colnames(pred.data))] <- predict(the.models[[3]], newdata = new.data1, type = "response", interval = "confidence", level = 0.95)
}else{
  new.data1 <- suppressWarnings(data.frame(x.levels, 
                                           df5[data.row, attr(the.models[[1]]$terms, "term.labels")[-c(1:3)]]))
  
  colnames(new.data1) <- c(c( t(do.call("rbind", lapply(the.models, function(x){attr(x$terms, "term.labels")[1:2]} ))) ), attr(the.models[[1]]$terms, "term.labels")[-c(1:3)] )
  
  pred.data[,grepl("\\.inc", colnames(pred.data))] <- predict(the.models[[1]], newdata = new.data1, type = "response", interval = "confidence", level = 0.95)
  pred.data[,grepl("\\.env", colnames(pred.data))] <- predict(the.models[[2]], newdata = new.data1, type = "response", interval = "confidence", level = 0.95)
  pred.data[,grepl("\\.edu", colnames(pred.data))] <- predict(the.models[[3]], newdata = new.data1, type = "response", interval = "confidence", level = 0.95)
}

return(pred.data)
}



boot.fun <- function(n.sims, the.models, x.levels, int = FALSE){
set.seed(1211)
pred.data <- t( replicate(n.sims, boot.fun.internal(the.models, x.levels, int) ) )
pred.data <- as.data.frame( apply(pred.data, 2, unlist) )
if(int == FALSE){
  pred.plot.data <- aggregate(pred.data[,grepl("fit\\.|lwr\\.|upr\\.", colnames(pred.data))], 
                              list(income = pred.data[,"income"],  
                                   envir = pred.data[,"envir"], 
                                   educ = pred.data[,"educ"]), 
                              mean)
}else{
  pred.plot.data <- aggregate(pred.data[,grepl("fit\\.|lwr\\.|upr\\.", colnames(pred.data))], 
                              list(income = pred.data[,"income"], income.neighb = pred.data[,"income.neighb"], 
                                   envir = pred.data[,"envir"], envir.neighb = pred.data[,"envir.neighb"], 
                                   educ = pred.data[,"educ"], educ.neighb = pred.data[,"educ.neighb"]), 
                              mean)
}
pred.plot.data[,grepl("fit\\.|lwr\\.|upr\\.", colnames(pred.data))] <- apply(pred.plot.data[,grepl("fit\\.|lwr\\.|upr\\.", colnames(pred.data))], 2, exp)

return(pred.plot.data)
}










